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Abstract — Mobile device localization in wireless sensor net- 
works is a challenging task. It has already been addressed 
when the WiFI propagation maps of the access points are 
modeled deterministically. However, this procedure does not 
take into account the environmental dynamics and also assumes 
an offline human training calibration. In this paper, the maps 
are made of an average indoor propagation model combined 
with a perturbation field which represents the influence of 
the environment. This perturbation field is embedded with a 
prior distribution. The device localization is dealt with using 
Sequential Monte Carlo methods and relies on the estimation 
of the propagation maps. This inference task is performed 
online, i.e. using the observations sequentially, with a recently 
proposed online Expectation Maximization based algorithm. The 
performance of the algorithm are illustrated through Monte 
Carlo experiments. 

Index Terms — Simultaneous localization and Mapping, Indoor 
localization. Received Signal Strength Indicator, WiFi, Signal 
Propagation. 



I. Introduction 

Wireless sensor networks |1| generally consist of a data 
acquisition network and a data distribution network, monitored 
and controlled by a management center. These networks have 
many applications such as environmental monitoring (|2|) or 
target tracking (lO, lU, Q). In this paper, we consider a WiFi 
communication network made up of a mobile device (such as 
a hand held mobile computer or a smartphone), a server and 
WiFi access points (APs). We are interested in the estimation 
of the localization of the mobile device in the environment 
using the signal strength of the surrounding APs. The mobile 
device collects the power of the signals and sends the data 
to the server which uses them to build an estimator of the 
device's position. The key step to provide such an estimator 
is to understand the behaviour of the WiFi signal strength for 
different positions in the environment. However, predicting 
the propagation of WiFi signals in an indoor environment is 
challenging since they are subject to many perturbations {e.g. 
shadowing, reflection...). 

Two main techniques exist to approximate the WiFi signal 
propagation map of each AP: the first ones use deterministic 
models based on the localization and characteristics of the 
surrounding APs as well as the localization of the obstacles 
involved in the environment, see for instance ||6|. Other 
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famous techniques are based on a previous hand made offline 
training phase in which a human operator performs a site 
survey by measuring the received signal strength indicator 
(RSSI) from different APs at some fixed sampled points, 
see |4l, (21. However, representing the indoor propagation 
map using a deterministic model is challenging since several 
obstacles cannot be taken into account. On the contrary, the 
site survey method allows to build an accurate estimation of 
the signal strength, but only for a finite number of points. 
Nevertheless, 1 8 1 provides a method to extend these measures 
to the entire map using Gaussian processes techniques. 

In this paper, we propose an estimation method that does not 
require any calibration procedure. The propagation maps are 
estimated online {i.e. without storing the observations) using 
the data sent by the mobile device. Any modification in the 
way the WiFi signals propagate inside the environment (due to 
new obstacles for instance) affects the data sent by the mobile 
device. Then, while these changes deteriorate the accuracy of 
localization systems using fixed estimators of the propagation 
maps, our system learns these changes by taking them into 
account in the construction of our map estimators. Thus, as 
illustrated in Section |V-B| the accuracy of our localization 



method improves with time instead of degrading. 

A semiparametric statistical model is used: the propagation 
maps are made of a parametric average indoor model in 
addition to a non parametric perturbation field. This model 
combines a prior knowledge on the signal propagation with 
random perturbations due to the obstacles. Based on the data 
collected by the mobile device, parameters and perturbation 
field estimators can be defined. We simultaneously provide 
an estimator for the device position. The procedure relies on 
an online Expectation-Maximization (EM) based algorithm for 
the estimation of the propagation maps and on particle filtering 
for the estimation of the device position. 

The structure of this paper is the following. Section In 
describes the model and defines the notations. Section [ill 
presents the online EM algorithm and Section llV] gives a 



general algorithm for online inference in our Simultaneous 
localization and Mapping (SLAM) problem. Section |V] illus- 
trates this algorithm with numerical experiments. 

II. Model and assumption 

Let {Xt}t>i be the cartesian coordinates of the mobile 
device in a two-dimensional compact space. This continuous 
environment is discretized into a finite grid map, denoted by 
C, for purposes of numerical computation. It is assumed that 
{Xt}t>i is a Markov chain taking values in C with initial 
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distribution v and Markov transition matrix given, for all 

{x,x') e by 



oc e 



(1) 



where a G W!^ depends on the average speed of the mobile and 
is assumed to be known and || • || denotes the usual euclidean 
norm in M?. Let B be the number of APs, |C| be the cardinality 
of C and be the B x \C\ matrix where Fj^^ is the j-th AP 
expected signal strength at position x. At each time step t, the 
mobile device measures and sends to the server the observation 
Yt taking values in R^. For all t > 0, the observation Yf is 
given by 

(2) 



where F^^j^^ = {Fj^Xtlj^i where {£t}t>o is a sequence 
of i.i.d Gaussian random vectors, independent from {Xt}t>i, 
with mean and covariance matrix S =^ cf^'^Ib (Ib is the 
identity matrix of size B x B). 

The position of the B APs are assumed to be known 
and denoted by {O^}^^. In order to take into account the 
perturbations in the signal propagation (due to the fact that 
radio waves are prone to shadowing, reflections and so on), 
we propose the following decomposition of F^: for all x G C 
and aU j G {!,••• ,5}, 



F^ =^/i* 



(3) 



For any B x 
A, 



J^l matrix A, we use the shorthand notation 
for the vector {Aj^x}xec- For any j G {1, • ' ' ^B}^ is 
the average indoor propagation and is such that for all x G C, 
/i^^ only depends on the distance between x and Oj. In the 
sequel, we use the so-called Friis transmission equation, see 
given by. 



★ def 



log \\x ■ 



■Oj\ 



(4) 



where c\ j and c^j are parameters depending on the environ- 
ment and log is the logarithm to the base e. 

Sj is an additive term due to random perturbations such as 
walls effects (a similar model of WiFi propagation maps using 
Gaussian processes can be found in 1 8 1). It is assumed that the 
parameters {6j}f^-^ are embedded with the prior distribution 
TT given, for any S G R^ '^l, by 



7r(^) cx exp 



where is assumed to be known and where, for any matrix 
A, A^ denotes the transpose of A. Figure ([T]) represents 
(sampled from tt) and the functions /i^ and Fj, defined on the 
grid C = {0, . . . , 30} x {0, . . . , 30}. The parameters used in 
this figure are Oj = (15, 15), and c^j, j and are given 
in Section |Vj their values were calibrated after a measurement 
campaign in an office environment. 

In the sequel, we write 0^ ^= (c^, c^^ S'^ ^ cr"^'^), where 



def 



def 



{c*}?, and S* {5*} 



For 

any x G C, the distribution of Yt conditionally to Xf = x has 
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Figure 1: Example of (sampled from tt) and spatial 
representations of the functions /i* and Fj = /ij + (5J (in 
dBm) 
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a density with respect to the Lebesgue measure on given, 
for all^ ='(^i,--- GR^,by 

B 



Therefore, {X^}^>o is the hidden process of a hidden 
Markov model observed through the process {Yt}t>{). The 
estimation of the mobile device's position Xt relies on the 
knowledge of the map . The observations {Yt}t>i are used 
to estimate simultaneously the mobile device's position and 0^ . 
This simultaneous localization and mapping problem may be 
seen as an instance of inference in hidden Markov models. 
For any positive integer n, any observation set . . . , y^), 
shortly denoted by yi^n and any parameter 6 = (ci, C2, (5, a^), 
the likelihood of the observations Lo{yi:n) is given by: 



Le{yi:n) =^ ^ iy{xi)g0{xi,yi)Y[q{xt-i^xt)ge{xt,yt) • 

Xi-r^eC'^ t = 2 

(5) 

Let n be a positive integer and Yi:n be a set of observations, we 
set as the estimator of 9*, the maximum a posteriori estimator 
defined as aTgma.Xgn~^£g{Yi:n), where: 



i0{Yl:n) = l0gie(yi:„) + log 7r(5) . 



(6) 



The next section provides a description of the EM algorithm 
and of online EM algorithms for the computation of max- 
imum likelihood estimators (without the penalty term). In 



Section IV we explain how such techniques can be used in 
our framework. 



III. Online EM 

The EM algorithm is a well-known iterative algorithm to 
perform maximum likelihood estimation in hidden Markov 
models |10|. Each iteration of this algorithm consists in a E- 
step where the expectation of the complete data log-likelihood 
(log of the joint distribution of the states and the observations) 
conditionally to the observations is computed; and a M-step, 
which updates the parameter estimate. 

Let Yi:n be a fixed set of observations and 6 be the current 
parameter estimate. 

i) The E-step consists in evaluating the conditional expec- 
tation 



1 



(J) 



where logp6i(Xi:^, Yi:^) is the complete data log- 
likelihood and [•lyi:^] is the conditional expectation 
given Yi-n when the parameter's value is 6. 
ii) The M-step updates the current value 6 taking the param- 
eter 6 maximizing (|7]). 

The model presented in Section |Il| belongs to the curved ex- 
ponential family: there exist functions S -.X^ xY ^ S dW^ , 
: e ^ R and : e ^ R'^ such that 

\ogq{x,x) ^\ogge{x ,y) = (j){0) + {S{x,x/ ,y),ilj{0)) , 



where (•, •) denotes the scalar product on R^. Moreover, there 
exists a continuous function : S ^ Q s.t. for any s e S, 

0{s) = argmax^^e {^W + {^.^{0))} • 

In this case the intermediate quantity defined by ^ can be 
written 



Q§{Yl:n'-, < 



m 

n 

n ^ — ^ 



t-1 



.XuYt) 



Yi, 



(8) 



Therefore, the E-steps amounts to computing only one con- 
ditional expectation Ylt=i S{Xt-i, Xt,Yt)\Yi:n] when 
the current parameter's value is 0. The M-step relies simply 
on the evaluation of at this conditional expectation. This 
two steps process is repeated till convergence. However, 
when the observations are obtained sequentially or when the 
E-step relies on a large data set, the EM algorithm might 
become impractical. Online variants of the EM algorithm 
have been proposed to obtain parameter estimates each time 
a new observation is available. In the case of independent 
and identically distributed (i.i.d.) observations, |11| proposed 
the first EM based online algorithm. The E-step amounts to 
computing intermediate quantities known as sufficient statistics 
(see below for a explicit definition) and [11 J proposed to 
replace these computations by a stochastic approximation step. 
When both the observations and the states take a finite number 
of values (resp. when the state-space is finite) an online EM- 
based algorithm was proposed by 1 12| (resp. by 1 131 ). These 
algorithms combine an online approximation of the filtering 
distributions of the hidden states and a stochastic approxima- 
tion step to compute an online approximation of the sufficient 
statistics. This has been extended to the case of general 
state-space models with Sequential Monte Carlo algorithms 
(see |[T4l . ifTSJ and 1 16|). More recently, 1 17 1 proposed a block 
online algorithm in which the parameter estimate is kept fixed 
on block of observations. The parameter's update then occurs 
at the end of each block. 

In this paper, we use the online variant of the EM algorithm 
introduced in |17| to perform the parameter estimation and 
to solve the localization problem presented above. This 
algorithm, called the Block Online EM (BOEM) algorithm 
relies on the ability to compute sequentially quantities of the 
form: 



— ''^^ S{Xt+n-l, Xt+n, Yt+n 



n:n-\-T 



Such quantities are called sufficient statistics. The BOEM 
algorithm uses a sequence of block-sizes {Tk}k>Q- Define 
To =^ and, for any A: > 0, =^ EzLi^i- Let k 
be a positive integer. Within each block of observations 
^Tfc+i:Tfc+i. the parameter's value 6^ is kept fixed and the 
sufficient statistic STfc,rfc+i (^/c) computed sequentially. The 
estimate 9^+1 is computed at the end of the block Yt^j^i-.Tlj^i 
through the evaluation of the function 6. 

Unlike the traditional use of the EM algorithm where 
the conditional expectations are computed using forward- 
backward techniques, QH, ITSi and ifTTll rely on recursive 
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computations of the conditional expectations. Indeed, let ^^,0 Algorithm 1 Bootstrap_filter_recursion (BFR) 
denotes the filtering distribution of Xn given the observations Require: ^^uj^ Yt, 0. 

Yi-n when the parameter's value is 0: 

VX G X , (t)n,e{x) = Pe {Xn = x\Yi:n) • 

Following lHa, lHa, defining for all x G X and all (9 G 6, 



we have, 



1 

— y s{Xt_i^ Xt^Yt 

n ^ — ^ 



t=l 



1 " 

-y^S{Xt-i,XuYt] 



^=1 



Yv.n 



Proposition 1 of |[T3ll illustrates the usefulness of this decom- 
position. 

Proposition 1 (of (TTSl). 
Initialisation: 

For all X and all 6 e Q, 



_ u{x)ge{x,Yi) 

E iy{x')go{x',Y,) 

x'ec 

= . 



Recursion: 

For all t >2 and all x e C, 



E (l)t-i,e{x')q{x'x)ge{x,Yt) 
x'ec 

E ^t-iA^')q{x'x")ge{x",Yt) 



^ I ^s{x',x,Yt 
x'ec 

Pt-i,e{x') !>. 



(9) 



(10) 



(l)t-i,e{x')q{x',x) 

E (l)t-i,e{x'')q{x'',x) ' 
x"ec 



Except in simple models (linear Gaussian models and finite 
state-space HMM), this algorithm requires forward computa- 
tions which are not available in closed form and which have to 
be approximated, e.g. using sequential Monte Carlo methods 
(see |[T4ll . ifTSll ). In this case, (/)^,6> is approximated by weighted 

^ L AT 

samples {^f,Sf}^i such that ^^,6>(^) = E Sf(5^p(x). In the 

sequel, {^f }^i will be referred to as the particle set at time 
step t. Plugging this approximation in ([TO]) yields: 

\s{iuet.yt) + {i-\)pU 



N 



p 
Pt 



(11) 



where is the approximation of pt evaluated at . At each 
time step, the new population of particles is built from the pre- 
vious population using Algorithm [T]ref erred to as the bootstrap 
filter, see e.g. |18|. The Bootstrap filter combines sequential 
importance sampling and sampling importance resampling 
steps to produce a set of random particles with associated 
importance weights. Implementations of such procedures are 
detailed in d, d, d, GB. 



1: for p = 1 to A/' do 

2: Draw / in 1 , . . . , A^ with probabilities proportional to 

3: Sample ^r^^(e/_i,-). 
4: Setc^f a^,(^f,r,)- 

5: end for 

6: return {il^l}^ 



This leads to the Algorithm |2] presented below. Algorithm [2] 
is the adaptation of the BOEM to our model. It recursively 
updates the parameter 9 at the end of each block. The BOEM 
proposed in 1 17] also introduced an averaged estimate based 
on a weighted mean of all the sufficient statistics computed 
in the past. It is proved in ifTVll that this averaged estimator 
has an optimal rate of convergence. Lines [211 to [25l of 



Algorithm [2] computes this averaged sufficient statistics and 
line [26] computes the sequence {Ok}k>Q of map estimates 
based on the averaged statistics. The BOEM algorithm is 
adapted by introducing a second particle system {^f ,cDf}^i. 
This additional particle system is generated using the averaged 
parameter estimate. As this estimate is supposed to be more 
accurate than the original estimator computed on each block, 
we use the second system of particles to build a better 
estimator of the device's position. At each time step, we then 
compute two estimators of the device's position, one for each 
particle system. Both of them are set as the particle with the 
greatest importance weight. LinefTS] performs the update of 
the sufficient statistics and line |19| the parameter's update at 
the end of the block. 

30| of Algorithm [2] we add a 



Finally, in lines 28 



to 



stabilization step (which is not in the original BOEM) which 
only consists in regularly replace the original map estimate 
by the averaged one. This step is needed to ensure the 
convergence as detailed in Section |V| 



IV. Application of the algorithms to the SLAM in 

WIRELESS NETWORKS 

In our framework, the objective is the maximisation of the 
penalized loglikelihood ([6]). This task can be performed using 
a similar technique as the one described in Section |IIl| since 
the additional penalty term only appears in the definition of 
the function 0. Define, for any {x^y) G C x and any 

Si{x) =^ {lx'{x)}^'^c , 
def 



S2,j{x,y) = {lx'{x)yj}x' 

( \ def 2 



The constant a being known, and since our model belongs 
to the curved exponential family, the penalized intermediate 
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Algorithm 2 BOEM_SLAM_indoor 



where, Ff =^ {F'-^}xec and 



Require: e^,^{Tk}k>i. {Yt}t>Q, N, . 

1: Set e = = 0^ . 

2: Sample {^olp^i {^^olp^i independently and uni- 
formly in C . 



3: Set 2; 



» — » 1 



— N 



for allp G {1,...,^"} . 



4: Set = for all p G {1, . . . , TV}, A: = 1, To = 0, 
Ti =ri. 

5: for alH > 1 do 



6: 
7: 



9: 
10: 

11: 

12 
13 
14 
15 

16 

17: 
18: 



19 
20 
21 

22: 
23: 
24: 



25: 
26: 

27: 
28: 
29 
30: 
31 
32: 



Selection and propagation step. 
Set {llQl}^^, = BFR ({ll_,M-i}ti^Yue] 

Set {i',S5r},'Li = BFR {{iUM-i}ti^ytA 

Position estimations. ^ 
Set p = argmax and = . 

Set p = argmax cDf and = . 

pG{l,...,Ar} 

Forward computation of the intermediate quantity. 
for p = 1 to X do 



Compute {pDp^i following (11). 
end for 

Map estimation. 
if t = Tk then 
Set 



Sk 



N 



= 0{Sk,n) . 

Set = for allp G {1,...,X} . 
if k = 1 then 

Set Sk = Sk • 
else 

Set 

5fe = (^Tk-iSk-i 

end if 

6> = 0{Sk,Tk) . 

Stabilization step. 
if /c = mod X5 then 

Set = 
end if 

/c = A; + 1 and Tfe = Tk-i 
end if 



-Tk. 



33: end for 



quantity can be written, up to an additive constant, as: 
1 ^ B 



^ {S3,,-2(S2,,,F,) + (Si,F/)} ^^^^ 



Si '=^' -E^ 



and, for all j G {1, . . . , 5}, 



1 



Yv.n 



t=l 



For any S = (Si, {S2,,}f=i, {S3,,}f=i) G [0,1] 



X (R^) , we denote by 6>(S,n) one of the 
parameter = (ci, C2, ^, cr^) maximizing the expression: 



^ {S3,j-2 (S2.j,Fj) + (Si,f,^)} 
i=i 



2a2 





def 




def 




def 


Wi,j 


def 




def 




def 


dj 


def 





log7r((5) - — logcr^ . 
n 2 

As mentioned in Section |ll| for any j G {1, . . . , 5}, Fj is 

def — 

written as Fj = jij + . In these experiments, the whole set 
of parameters (c^, C2 and a^'^) and the unknown perturbation 
Gaussian fields {(^J}^^ are estimated using Algorithm [2] For 

all j G {1, . . . , B}, we write Dj =^ {log \\x - Oj \\}xec and 

2 

diag(Si) + ^S-i , 

diag(Si) [/-Mo-jdiag(Si)] , 
/ - diag(Si)Mo-i , 

Thus, 6 = (9(S, n) is given, by ^ = (ci, C2, (5, a^) where, for 
all jG{l,...,5}, 



= d-'[Ws,jl^-W2,jDj]M,,,S2,j. 
= dj^ [-W2,jt^ + Wi^jDj] M2S2,, , 
Sj = Moj [S2J - diag(Si)(ci,_^-l + C2jDj] 



and 



F, 



a = 



Ci .-1 + C2 .-I^.- 



1 ^ 



{F;diag(Si)F, -2Si^,^.F, +S3,,} . 



V. Experiments 



A. Simulated data 



2a2 



In this section, the performance of the proposed BOEM 
algorithm is illustrated with simulated data. All experiments 
are performed on the grid C = {0, . . . , 30} x {0, . . . , 30}. 
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We use B — 17 APs, each AP being modelled by the same 
coefficients and C2, see (|4]), 



VjG{l,...,5}, cl 



-26 and c 



-17.5 . 



For all j G {1, . . . , 5}, is a Gaussian co variance function 

defined by T>j{x^x') = vi ^ exp(— |x — x'\'^/v2) with 'Ui = 
10 and V2 = 18. The variance of the observation noise is 
cr"^'^ = 25. The variance of the transition kernel defined in ([T]) 
is chosen such that a = 6. 

All runs are started with the same initial estimates 0^ = 

(c?,cO,(^o,aO'2) where, = 

VjG{l,...,5}, c?.. = -10 



^0 



-30 and a 



,0,2 



30 



The number of particles N = 25 is kept fixed and the initial 
position of each particle is chosen randomly and uniformly 
in C. For each map Fj, the estimation error is set as the 
normalized Li error, such that the distance of a given map Fj 
from the true map is 



J,x I 



and the error displayed is the mean over all maps: 

B 



.def 1 

■ =1 



J ' 



The block sizes are given by 



lOk + 500 . 



On each block, the localization error is set as the 0.8-quantile 
of the distance between the true localization and the estimated 
position. Figure |2] displays the error on the estimation of the 
maps and on the localization when the stabilization step in 
Algorithm [2] is omitted (lines 29 to 31). This case corresponds 
to the BOEM algorithm. The localization part is dealt with 
using two different procedures. 

• Nonaveraged estimate: the estimate is given with the 
original particle system (see line [TO] of Algorithm [2]). 

• Averaged estimate: the estimate is given with a second 
particle system run with the average estimation of the 
map (see line 11 of Algorithm [2]). 

In order to give fair results, the optimal estimate is shown, i.e. 
the estimated position given with a particle system run with 
the true maps F'- , j G {1, . . . , B}. 

As shown in Figure |2a| the estimated position does not 
converge as the number of blocks (i.e. as the number of 
estimations) increases. After 50 blocks (about 40000 ob- 
servations) the position, which is badly estimated, does not 
provide good map estimates which increases the error on the 
averaged map estimate. Figure |2b] displays the error on the 
map estimate. It is clear that both the estimate and its averaged 
version do not converge. This convergence problem of the 
BOEM algorithm can be due to the curse of dimensionality 
that can occur when the number of parameters to estimate is 
high. Moreover, the higher the parameter space dimension is, 
the more likely EM based algorithms are prone to converge 
towards local minima (see (HI). To overcome this difficulty. 
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(a) 0.8-quantile of the distance between the true locaHzation and the estimated 
position. The locaHzation error is given with the nonaveraged estimate (dotted 
line), the averaged estimate (dashed line) and the optimal estimate (bold line). 
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(b) Mean Li error on the map estimate with the initial estimate (dotted line) 
and the averaged estimate (dashed line). 

Figure 2: Errors on the map estimation and localization 
processes with the original algorithm. 



we propose to use the good behaviour of the averaged map 
estimate during the first 50 blocks. The map estimate is 

|30l 



regularly replaced by its averaged version, see lines 28 



to 



of Algorithm [2] This will prevent the map estimate from 
diverging and thus, this will reduce the error on the estimated 
position. In Figure |3j this stabilization process is performed 
each time A/'^ = 5 blocks have been used. As shown by 



Figure [3a| and Figure [3b] this greatly improves the performance 
of the estimation of both the maps and the localization. 
Hence, the proposed algorithm is based on this stabilization 
procedure and uses the averaged position estimate to perform 
the localization part. 

Figures |4] and [5] illustrate the performance of the algorithm 
for the localization and for the estimation of the maps over 
50 independent Monte Carlo runs. In Figure |4] the optimal 
localization error (i.e. when the maps are known) is also 
displayed. The convergence of the localization error to the 
optimal error is almost reached after 100 blocks (about 100000 
observations). Similarly, the error for the estimation of the 
maps given by the averaged algorithm goes on decreasing after 
100 blocks (the decrease is slower after 75 blocks). 

B. True data 

In this section, the behaviour of our SLAM algorithm is 
illustrated in a real situation. 10 access points are set up in 
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(a) 0.8-quantile of the distance between the true locaHzation and the estimated 
position with the stabilization process. The locaHzation error is given with the 
nonaveraged estimate (dotted line), the averaged estimate (dashed line) and the 
optimal estimate (bold line). 
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Figure 5: Boxplots of the mean Li error on the map estimate 
with the stabihzed algorithm and the averaged estimate as a 
function of the number of blocks. 



(b) Mean Li error on the map estimate with the initial estimate (dotted line) 
and the averaged estimate (dashed line) with the stabilization process. 

Figure 3: Errors on the map estimation and localization 
procedure with the stabilized algorithm. 




Figure 6: Map of the indoor environment used for the test with 
the position the access points (red circles) and their associated 
identification numbers. 




Numbei^§f blocks 



Figure 4: Boxplots of the localization error given by the 
stabilized algorithm with the averaged estimate (left) and the 
optimal estimate (right) as a function of the number of blocks. 



an office environment (Figure [6] represents a map of this 
environment as well as the position of the access points). 
The map is discretized using a grid C C [0,30] x [0,30]. 
The variance a"^'^ is assumed to be known and its value 
(cr^'^ = 2bdBrin?) is calibrated using a measurement campaign 
at a fixed position. Around T = 20000 measures of the RSSI 
have been made on the map using a WiFi device. Algorithm [2] 
produces position estimates but we do not have a direct access 
to the real position and thus cannot observe the localization 
error. To overcome this difficulty, a test data sample is built by 
producing measures along 12 paths in the environment such 
that, for each measure, the associated position is registered. 
The test data sample is made of Ttest = HOO measures: 
jjj^test y^testj^st^ The tcst data sample is used to compare 
the localization accuracy provided by different values of the 
parameter F. Note a major difference between the model 
given in Section |Il| and the real data situation. For any 
measure Y sent by the device, only several APs are represented 
in Y. Therefore, the maps Fj, j G are not 

estimated simultaneously as, for any time step t, two APs 
might appear a different number of times in Yi:t. We thus 
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Number of updates 

Figure 7: 0.8-quantile of the distance (in meter) between the 
true locaHzation and the averaged estimate obtained with the 
stabiHzed algorithm. The locaHzation error is computed on 
the test data sample each time one of the estimated map is 
updated. 



slightly modify Algorithm |2] by introducing specific blocks 
and measure counters relatively to each AP. Each time the 
value of Fj for any j G is updated using a 

block of the 20000 measures, we submit the new estimator 
F to the test data sample: Algorithm [2] is run on the test 
data sample {Y^^^^}J^\. Only the averaged particle system 
is computed and no parameter update is performed. We can 
then compute the localization error relatively to the test data 
sample as the 0.8-quantile of the error between {Xl^^^}J^\ and 
the averaged position estimate. Figure [7] displays the results of 
this experiment by representing the 0.8-quantile of the error 
as a function of the number of updates. The numbers on 
the graph in Figure [7] indicate which AP were updated for 
each update. The initial map estimates are given, for any 



jG{l,...,10} by 



-26 and c^ j 



-17.5 and So = 0. 



Despite the relatively small test sample size. Figure |7] shows 
that the localization error seems to adopt the same behaviour 
as the localization error for the simulated data. The parameter 
Fj was updated a maximum of 7 times (for AP j = 10 for 
instance) and a minimum of 2 times (for AP j = 3). Figure 
[s] represents the final estimate of the propagation maps Fj , 
j e {1, . . . , 10}. For some map estimates (see for instance 
the access points 1, 4 and 7), the signal strength drops when 
passing walls while the walls are responsible for a part of the 
indoor waves propagation disturbances. 

VI. Conclusion 

In this paper we propose a stabilized version of the BOEM 
algorithm to estimate the signal propagation maps needed in 
any WiFi based localization system. The main difference 
with the existing solutions is that these propagation maps are 
estimated using the data sent by the mobile device originally 
used for localization purposes. On the contrary, the existing 
WiFi based localization systems establish these propagation 
maps either in a deterministic way or by running a previous 
hand made survey. In case of environmental modifications, 
the propagation maps are thereby changed. Our technique can 
easily be adapted to these changes by regularly reinitializing 
the sufficient statistics while hand made survey based systems 
can not take into account these modifications without renewing 
the survey. However, further tests are needed to evaluate the 



Ik. 





I, 



i: 




Figure 8: Graphical representation of the final propagation 
maps estimations: (in dBm). 



accuracy provided by our method and to compare it with 
other methods. Many elements should be analyzed such as 
the number and the position of the access points, the size of 
the environment or the materials constituting the obstacles in 
the environment. 
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